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ABSTRACT 


This two-part report is concerned with the development of a general framework 
for the implicit time-stepping integrators for the flow and evolution equations in 
generalized viscoplastic models. The primary goal is to present a complete theoretical 
formulation, and to address in detail the algorithmic and numerical analysis aspects 
involved in its finite element implementation, as well as to critically assess the numerical 
performance of the developed schemes in a comprehensive set of test cases. On the 
theoretical side, the general framework is developed on the basis of the unconditionally- 
stable, backward-Euler difference scheme as a starting point. Its mathematical structure is 
of sufficient generality to allow a unified treatment of different classes of viscoplastic 
models with internal variables. In particular, two specific models of this type, which are 
representatives of the present state-of-art in metal viscoplasticity, are considered in 
applications reported here; i.e., fully associative (GVIPS) and non-associative (NAV) 
models. The matrix forms developed for both these models are directly applicable for both 
initially isotropic and anisotropic materials, in general (three-dimensional) situations as 
well as subspace applications (i.e., plane stress/strain, axisymmetric, generalized plane 
stress in shells). On the computational side, issues related to efficiency and robustness are 
emphasized in developing the (local) iterative algorithm. In particular, closed-form 
expressions for residual vectors and (consistent) material tangent stiffness arrays are given 
explicitly for both GVIPS and NAV models, with their maximum sizes “optimized” to 
depend only on the number of independent stress components (but independent of the 
number of viscoplastic internal state parameters). Significant robustness of the local 
iterative solution is provided by complementing the basic Newton-Raphson scheme with a 
line-search strategy for convergence. In the present first part of the report, we focus on 
the theoretical developments, and discussions of the results of numerical-performance 
studies using the integration schemes for GVIPS and NAV models. 
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Robust Integration Schemes for Generalized Viscoplasticity 
with Internal-State Variables; Parti Theoretical 
Developments and Applications 

1. Introduction and Literature Review 

Dealing with the general topic of computational inelasticity, the scope of the 
work in this report is focused on the development of algorithms for the integration of rate 
dependent constitutive equations. In recent years, this subject has attracted considerable 
attention and will form the basis for the present work. In particular, recent literature has 
clearly emphasized the use of implicit integration schemes in view of their robustness; i.e., 
their superior stability and convergence properties for viscoplasticity, as well as for the 
limiting case of inviscid plasticity. 

In regards to the rate dependent constitutive models, a number of viscoplastic models 
were developed to represent the material behavior of metals and composites, etc. All 
viscoplastic models can be classified into two classes, the first is the fully associative 
viscoplastic model, or GVIPS ( Generalized Viscoplastic with Potential Structure ), the second 
is the NAV ( Non-associative Viscoplastic models ). GVEPS possess both the thermodynamic 
potential ( Gibb’s function ) and the dissipation function ( Q form ), while NAV generally refer 
to those models that have partially (e.g., f2 form only) or totally incomplete potential form. In 
these latter NAV models, thermodynamics requirements are simply reduced to the fundamental 
dissipation inequality only. For these two different classes of models, a general computational 
framework suitable for implementation of both is needed, and this constitutives the major part 
in the developments reported here. 

1.1 Integration Schemes 

Computational algorithms for the integration of constitutive relations play a central role 
in the inelastic finite element analysis of engineering structures. For example, because of their 
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simplicity, the classical von Mises (or Jr-type) theories of inviscid plasticity [28], and 
generalizations for its time-dependent viscoplastic counterpart [37, 58] are extensively utilized 
for isotropic metals. Correspondingly, a major research effort has been devoted over the years 
to the development and critical assessment of stress updating or integration schemes for the 
rate equations in these material models [2, 8, 12-14, 16, 22, 29, 31, 33-35, 40, 42, 44-45, 48- 
49, 52, 54, 57, 63, 67-70, 74, 76-78, 80], As a result, several powerful techniques are 
currently available; e.g. refer to the extensive reviews and evaluations in [51, 64] for plasticity 
and [2, 30, 35, 65, 78] for viscoplasticity. 

In particular, note that, in early applications the explicit integration schemes (i.e., 
forward Euler method) were predominate because of their ease in implementation, and 
because they do not require evaluating and inverting a Jacobian matrix. However, explicit 
integrators may not be efficient. That is, too many iteration steps may be required and 
convergence stability can not be guaranteed [49, 54, 80], More recent work has clearly 
emphasized the use of implicit integration methods [12-14, 22, 29, 31, 33, 40, 42, 45, 48, 52, 
57, 63, 65, 68-69, 74, 76-77] in view of their superior stability and convergence properties [30, 
51], Several alternative approaches have been used in the derivation and subsequent mathe- 
matical analyses of these latter methods (collectively known as return mapping algorithms in 
recent computational-plasticity literature). For instance, these include convex analysis tech- 
niques for variational inequalities of plasticity [23, 32, 46, 50], mathematical programming or 
holonomic methods for incrementally external paths [11, 18, 41, 43], the "substructure" 
analogy [59], asymptotic expansions for integral-equations in viscoplasticity [17, 25, 75], as 
well as the more conventional finite differencing schemes [26] and their generalizations; i.e., in 
the form of semi-implicit or forward gradient approach [e.g. 2, 57, 78], and also implicit 
trapezoidal or midpoint rules [30, 51, 64], Walker’s asymptotic integration scheme is of 
more recent origin, and exhibits a number of unique advantages which are marked 
improvement in accuracy and stability over existing integration algorithms for the case of 
isotropic materials subjected to uniaxial loads [17, 75], In these references, the asymptotic 
integrator is described for dissipative linear two dimensional systems of ordinary 
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differential equations, and is shown to yield the exact solution for a one dimensional linear 
equation. However, recent work shows that for higher dimensional systems of nonlinear 
equations, such as Freed’s viscoplastic model, accuracy and stability cannot be guaranteed 
in general [27], From the standpoint of practical applications, the one-step, fully implicit, 
backward Euler scheme is presently one of the most widely used integrators [14, 29-31, 33, 40, 
48, 51-52, 63-64, 68-69, 76-77]. For the computationally intensive viscoplastic applications, 
found in typical finite element analysis, implicit backward Euler integration methods have 
become the proven standard for the numerical integration of the viscoplastic rate 
equations [27, 65], 

1.2 Viscoplastic Models 

During the past decade, much progress has been made in the development of 
viscoplastic constitutive equations. Several viscoplastic models have been proposed and 
developed to treat the complex time dependent viscoplastic behavior of metals, alloys and 
composites at high temperature [15, 24, 60, 62]. The deformation behavior of materials at 
high temperature involves energy dissipation and material stiffness variations due to 
physical changes in the material’s microstructure. Consequently, thermodynamic 
arguments have often been utilized as a foundation on which phenomenological 
constitutive laws may be formulated. The general form for a fully-associative potential- 
based viscoplastic formulation in terms of the Gibb’s thermodynamic potential has been 
put forth in [3-5], The complete potential-based class of inelastic constitutive models 
exhibit a number of unique advantages from both a theoretical and a computational 
standpoint, for example, the symmetry of the resulting consistent tangent stiffness matrix, 
and possesses a form which is convenient for further development of new deformation and 
damage models. These advantages serve to motivate further development and use of this 
class of models. 
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Another class of constitutive models are those that are nonassociative and are 
currently used by industry, e.g., the unified viscoplastic theory of Robinson and Duffy 
[60,62] for metal matrix composites at high temperature, together with its extension for 
coupled creep-damage analysis in [61]; and Freed’s viscoplastic model with creep and 
plasticity bounds to analyze the response of polycrystalline metals at high temperature 
[24], Recent work has demonstrated that the above models may be modified to restore 
the complete potential structure. With regards to the computational algorithms, research 
work for the more complex constitutive theories is still needed. In this regard, the study of 
the important issues involved in such an undertaking constitutes our main objective in the 
present report. 


2. Scope, Objectives and Contribution 

2.1 Scope 

From the standpoint of practical applications of interest (i.e., large-scale numerical 
simulations), and in view of the above discussion, the focus of the present work is on the 
stress updating schemes of the one-step, fully implicit, backward Euler difference scheme. 
This is the simplest in its class, and is certainly the most widely used at present. Based on 
recent experience and results, this scheme was found to be efficient and stable when 
effectively utilized in a class of isotropic and anisotropic coupled viscoplastic-damage 
models. 

With regards to the constitutive models employed, two general classes of unified 
viscoplastic theories are considered: 

(a) The nonassociative models ( NAV ) 

(b) The fully associative, or complete potential-based, models ( GVIPS ) 
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Both of these models are sufficiently general for the anticipated applications; i.e., 
incorporating such important response characteristics as nonlinear (isotropic and/or 
kinematic) hardening, various recovery mechanisms, etc., and both are amenable to 
extensions for nonisothermal conditions, as well as anisotropy (i.e., generalized J 2 format). 
The majority of viscoplasticity models currently in use (at NASA and elsewhere) belong to 
class (a). Several developments in class (b) are currently underway as these models are 
very appealing in view of their theoretically sound basis from the thermodynamics 
standpoint, and the symmetries exhibited in the corresponding integrated forms of their 
flow and evolution equations. This report will discuss the present implementation which 
has been developed in terms of a framework which is applicable to both GVIPS and 
nonassociative models to show the robustness of the backward Euler integration 
algorithm. As an example of a nonassociative model, Freed’s model [24] was chosen for 
implementation into the proposed integration algorithm. For the generalized form of 
GVIPS, the decoupled potential framework is assumed and the specific GVIPS model is 
put forth. 

For conciseness, the discussion is limited to a case involving small deformations ( 
in which the initial state is assumed to be stress free throughout ), an initially isotropic 
material, and isothermal conditions. Note that the implementation framework is also valid 
for an anisotropic material. 

2.2 Objectives 

For each of the viscoplastic models considered, and the associated implicit 
backward Euler integration method, the following tasks are to be completed: 

(1) Detailed study of the underlying mathematical structure of the viscoplastic equations, 
and the corresponding integrated fields of stress and internal state variables. 

(2) Development and implementation of the stress-updating algorithm, and associated 
line search and subincrementing schemes. 
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(3) Testing of the convergence, stability, and accuracy properties of the algorithms. 

(4) Documentation of results, guidelines, and recommendations for effective utilization 
of the schemes developed. 

Our primary objective in item (1) is to identify the pertinent matrix operators 
needed in the algorithmic development; i.e., concise forms of the integrated field 
equations, and convergence limits for the subincrementing control strategy in item (2). 
The implementation in item (2) will encompass several nonlinear solution procedures; i.e., 
both the initial (constant) and tangent (variable) stiffness formats. The latter utilizes 
appropriate forms of the consistently-derived material tangent stiffness; i.e., 
model/integrator-dependent. 

The implementation and coding in (2) will be written in a “modular”, “stand- 
alone”, format with full documentation of the associated material-model subroutines. This 
will facilitate a straightforward utilization of the developed algorithms in other nonlinear 
finite element analysis codes. 

In conjunction with (3), an extensive number of discriminating test problems will 
be utilized. The objective is to enable firm conclusions to be made regarding the “relative” 
merits and/or limitations of various schemes. 

2.3 Contribution 

Based on the present work, some achievements are summarized as follows: 

(1) Original equations of Freed’s model are recast into the matrix format described in [65] 
to facilitate the model’s implementation into the newly developed frame-work 
discussed in Part II [38] of this work. 

(2) Freed’s nonassociative model is implemented successfully in finite element analysis 
for the generalized plane stress case. 

(3) The algorithm described in the report is proven to be accurate and efficient. 
Convergence and stability are guaranteed. 
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(4) Line search technique is used successfully to achieve robustness of the integration 
scheme. 


3. Mathematical Formulations of Viscoplastic Models 

3.1 Introduction 

In this section, mathematical formulations of the complete potential-based 
viscoplastic structure is discussed. Based on recent work [3, 5], the general form for 
GVIPS should have two key ingredients; 1). the thermodynamic potential, or Gibb’s 
function ( or conjugate free energy, Helmholtz function ) defining the equations of state, 
and 2). the viscoplastic dissipation function ( G form ) for the ensuing flow and evolution 
laws of inelastic strain and internal state parameters. The previous forms typically 
employed by many researchers [60, 62, 24] in discussing the general structure of the 
“thermodynamically based” constitutive equations are /2-form ( flow or dissipation 
potential ). These forms are introduced solely for convenience in satisfying the 
thermodynamic admissibility restriction of the engineering materials based on simple 
properties of non-negativeness and convexity of these functions. From a strict 
mathematical standpoint, the “thermodynamic-admissibility” restriction associated with the 
dissipative mechanisms reduce to the well-known Clausius-Duhem or dissipative 
inequality (/2-form) [19, 39], However, such /2-forms do not presuppose or 
automatically imply the existence of the total ( integrated ) forms of the associated 
thermodynamic potentials, e.g., the ( Helmholtz ) free energy, or Gibb’s function. 
Together, the /2-form and the a priori assumption for the existence of Gibbs and 
Helmholtz functions will lead to a complete potential-based formulation. For distinction, 
the original model, i.e., using only the /2 -forms, is referred to as the incomplete potential 
formulation or the so-called nonassociative model. 
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The complete potential-based class of inelastic constitutive equations possess a 
number of distinct and important attributes from both a theoretical and a computational 
standpoint. First, they constitute the cornerstone of numerous regularity properties and 
bounding ( or limit ) theorems in plasticity and viscoplasticity. Second, they result in a 
sufficiently general variational structure, whose properties can be exploited to derive a 
number of useful material conservation laws. Third, on the numerical side, the discrete 
form of this same variational structure is of great advantage in the development of efficient 
algorithms for finite element implementation, for example, symmetry-preserving material 
tangent stiffness operators are easily obtained in implicit solution schemes (examples are 
given later in Part II of this report). Finally, this complete potential-based framework 
conveniently lends itself to intelligent application of symbolic manipulation systems that 
facilitate the construction, implementation, and analysis of new deformation and damage 
models [3], Several of the currently used forms [60, 62, 24] do not conform to the general 
potential framework, but it is still possible in some cases, by modifications of the 
employed forms, to restore the complete potential structure. 

3.2 Generalized Viscoplasticity With Potential Structure ( GVIPS ) 

The discussion is limited to the case of small deformations, in which the initial state is 
assumed to be stress-free throughout. For generality, the initial state of material isotropy 
or anisotropy is left unspecified in this section; i.e., except for the explicit appearance of 
the corresponding material-directionality tensors in the arguments list for the potential 
functions employed, all derivations remain generally applicable. 

Consider the following general form for the Gibb’s potential: 

'F = 'F(cr*,<4,£>‘') (3.1) 

where a b T denotes the T th-order tensorial internal state variables (where b=l,...,nb ; T= 
l,...,nx), a b T can be a second-order tensor (back stress) or higher order tensors, b denotes 
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the number of tensorial internal state variables; D d is the scalar internal state variable 
( where d=l,...,nd ) and <j tJ the cauchy stress tensor. This is a typical general form for the 

Gibb’s potential, in which as many internal variables as desired may be specified in eq. 
(3.1). The potential may be decoupled into elastic and inelastic parts: 


¥ = '¥‘(cr ij ) + x ¥ i (a b T ,D d ) 

(3.2) 

or 


T - 4"( ct# ) + £X(<4) + i'Pi(Z)') 

6=1 d= 1 

(3.3) 

It follows from eq. (3.3) then that: 


, 

y 

(3.4) 

II 

(3.5) 

II 

T5 Q 

(3.6) 


where e tj is the total strain and^£ is inelastic strain tensor, respectively. 

These above relations are defined as the equations of state, and a tj , a b , D d are the 
“force-like” thermodynamic state variables while £y, y b T , y d D are the conjugate 
“displacement-like” variables ( affinities ). 

Now, assume the existence of a dissipation potential of the form. 


n = 0(0,, a\, D d ) 


(3.7) 
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Which when written in its decoupled form [3] gives, 

a = i2,(or) + + £ «?(!>'') 


(3.8) 


In order to provide a framework for the definition of the inelastic strain rate (e£) and 
ensure its thermodynamic admissibility, the Clausius-Duhem or dissipation inequality is 
required. That is. 


a(q„ a‘ T , D 1 ) - < t ,«{ + * 0 

b 


(3.9) 


The rate of change of the affinities may then be expressed in terms of their corresponding 
internal state variables. That is, the flow or kinetic law becomes: 


H = 


dQ 

da, 


. b dQ 


7 d = 


an 

3D‘ 


(3.10) 


(3.11) 


(3.12) 


From eqs. (3.5) and (3.6), 



1 

d 2 *P 

Yt dt\ 

M\ 

da b T da b T 

d 

[ d*F~ 

d 1 ^ 

Yd dt 

L dD d _ 

cD d cD d 


■D d 


(3.13) 


(3.14) 


then, from eqs. (3.11) (3.12) (3.13) (3.14), 


11 


&T [ * b * b 

da T da T 


D d = [ 


cD d cD d 


! dQ 
da\da b T da\ 

(3.15a) 

! en 

'~<3D d <%) d cD d 

(3.15b) 


which defines the evolution of the internal state. Eqs. (3.10) and (3.15) represent the flow 
and evolution laws, respectively. Both the state equations (3.4), (3.5), (3.6), and 
evolutionary laws must be satisfied for the general forms of GVTPS which has an assumed 
dissipation function Q = /2(cr, a b , D d ) and Gibb’s potential wherein both potentials are 

directly linked through the internal state variables a b T and D d (see eq. 3.15). 

From eq. (3.15), the compliance operators are defined as follows: 
r _ t ctY 

L *-~da b T da b T 5 dD d dD d ( ' 

they relate the “force-like” state variables ( a b ,D d ) to the “displacement-like” variables, 
provide a link between the assumed Gibb’s and complementary dissipation potential, and 
ensures a number of desirable numerical features such as the symmetry of the resulting 
consistent tangent stiffness matrix. The above discussion provides the general 
thermodynamic framework in which the flow and evolution laws are associative. 


3.3 Specific Form for GVTPS Model 

Next, we use a specified form to constitute a GVIPS model. The emphasis is placed 
here on a careful examination of the mathematical structure of these governing equations. 
The transversely-isotropic viscoplastic material, together with the counterpart anisotropic 
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elastic-response component, for a solid reinforced by a single family of fibers as in [71- 
72], is assumed. The development and application of fairly comprehensive models for 
anisotropic materials may be referred to [1, 6, 9-10, 21, 47, 53, 55-56, 60-62, 73], 

3.3.1 Elastic Response 

The general viscoplastic models consist of an elastic response and a viscoplastic 
response. The hyperelastic response is assumed to be linear . ; i.e., a quadratic strain-energy- 
density function. If, is assumed to exist such that the Cauchy (true) stress components a are 
given by 


cr„ — _ _ — C iM e 


ds e . 

V 


ijkl^kl 


where a superscript "e" signifies an elastic component; and 


(3.17) 


e = e e + e p . 

'J ‘1 >] 


(3.18) 


in accordance with the basic hypothesis of the additive decomposition of total strain tensor, e 9 , 
into its elastic, £*, and inelastic strain components, e p . In the above, Cfju denotes the fourth- 
order , symmetric, elastic moduli tensor; i.e. C'jki = Cuij. 


3.3.2 The Viscoplastic Response 

Now, specify that T=2, nb=l, nd=0, so only one second-order internal variable, which is 
the back stress (internal), is specified. An outline of the basic equations is therefore given here. 
The starting point is a postulate for the existence of a plastic potential function, Q, in terms of 
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the invariants of the stress tensor, a, the internal variable (kinematic-hardening or back stress 
tensor), a, such that 


Q = ^(cr^a^Vy ) = C2, (a, ,a,,V ¥ )+ G 2 (a^V # ) 

(3.19a) 

n=n,(F)+n 2 (G) 

(3.19b) 

an 

£ P . = 

9 da 

IJ 

(3.20a) 

d£l 

/o 'da,., 

(3.20b) 


where the parentheses are used to indicate functional dependence on the arguments inside (•); 

and the overdot signifies a rate or time derivative. 

The F and G are quadratic functions in terms of the invariants of the deviatoric 
components of the effective stress, (CTij-otij), and the back stress, Oij, tensors, respectively. 


F = ±(<7 s -a,)Me,(<r u - au )/X?.l 

(3.21a) 

G = \{a,)M m (a u )l K] 

(3.21b) 


in which the symmetric fourth-order tensor Myki is a function of fiber orientation tensor V, ( 
material-directionality ) and is defined as: 

(3.22) 


where 
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II 

M 1 
M • 

4(c J -l) 
^ 4<a 2 - 1 

(3.23a) 

r|=K)/BQ ; 

o=Yi/Y t 

(3.23b) 

and Pijki, Qijti and Rjju are defined as follows: 



Pylcl = ^ijkl — > I.jkJ 

= |(8,5 i , + 5 J 5 lk ) 

(3.24a) 


Q* = Ifv.Sj, + V.8* +V„5„+ V/, ) - 2V„V„ (3.24b) 

Rp =3V,V U -(V # S U +V u 5 ij )+^8j j S u (3.24c) 

Vj-ViVj (3.24d) 

where v; is the unit vector defining the material fibre direction. 

In the above, 0 < £ < 1 and 0 < C, < 1 are material-strength ratios, and the constants K» ( 
or K) ) and Y t ( or Yi ) are threshold strengths of the composite material in transverse ( or 
longitudinal ) shear and tension/compression, respectively. For the isotropic case, £ = £ = 0 and 
Mjjid reduces to P.ju, i.e., simply the classical von MIses-type forms in functions F and G. 

Because of nb=l, n,r=0, for the present purely mechanical theory, Eq.(3.1) comes to the 
following decoupled form of a (Gibbs) complementary free-energy density function in terms of 
stress and internal state parameters [e.g. 39]; i.e., with superscripts "e" and "i" indicating elastic 
and inelastic parts, 

^'(s) ( 3 - 25a ) 


we define vj/ e and V by the following forms: 
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'?'(<*«) = \<ri, C w 


¥ *( a ) = Q 1+ fi 

' v) (1 + P)H 


(3.25b) 


and the (positive) material constants: H and K, (stress units), and P (dimensionless), as well as 
function G = G(ajj) were defined previously in Eqs. (3.21b). Together with the strain vector 

decomposition ofEq. (3.18), the above Eqs. (3.25) define the equations of state: 


r r p — — C* 1 rr 

av i 

(3.26a) 

iJ iJ da m u 9 

UO ij 

Y * ax, h K « 

h — H/G** * 7Cy — Myki ocu 


(3.26b) 


the dissipation function D (<Jij,<Xjj) is now introduced naturally as follows: 

n =a, e ,-\ \i > o (3.27a) 


or, substituting for vji in terms of a, and a, from (3.25a), one gets from (3.27a) 


dn 

da 

y 


= e.v 


BQ. 

ax. 


= -r, = 


aV 


fa-ifaki j 


a, 


(3.27b) 


The rate form presently connecting a, and y, through the Hessian matrix of v/ ( hardening- 
moduli tensor) ,the time derivatives of the second of Eqs. (3.27), yields 


h 9 

<x„ = h Z”L + a,, a 

9 I 1 h(l+ 2$) J 


Yi/ [Aj«] Y Id 


(3.28a) 


where L ijU is called compliance operator, and h’ is a scaled derivative, 


h’ = 


1 dh 
k t 2 dG 


(3.28b) 
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and we defined Z” u = for the “generalized inverse” of M\ see details for different 

conditions (three-dimentional, plane-strain/axisymmetric, generalized plane stress, etc) in Part 
II [38], Now, the present model exhibits a potential structure, i.e., the elastic part is formulated 


flow law for the inelastic strain rate and the nonlinear kinematic hardening-recovery type 
evolution equation for the internal state variable. 

3.3.3 The Viscoplastic Response: Flow and Evolution Laws 

From Eqs. (3.19-3.21) and (3.25-3.28), the final expressions for the inelastic flow and 
evolution equations are given as follows, for the specific selection of power functional forms 
for Qi (F) and 0.2 (G): 


where/ h, and y are the flow, hardening, and recovery functions, respectively. The latter are 
given in terms of their respective arguments as follows: 


in which n > 1, p, H, P, R, and m > P + 1/2 are (positive) material constants, and K, was used 


in terms of a strain-energy density function, whereas a dissipation potential is utilized in the 



(3.29a) 



n ij ~ kl 


(3.29b) 



(3.29c) 


f(F) = F n /(2pK t >/F+I) 


(3.30a) 



(3.30b) 


previously in Eqs. (3.21) to (3.25). 
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Several regions of discontinuities in the state (cjjj, 0 Cij) space were introduced to account 
for stress-reversal/ cyclic-loading/ dynamic-recovery effects. This simply amounts to using 
different functional forms for / h, and y in different regions; i.e., 


/ = / if F > 0 

(3.31a) 

/ = 0 otherwise 

(3.31b) 

G = G if G > Go and dev^ - > 0 

(3.32a) 

G = Go otherwise 

(3.32b) 


where the "small" constant (cut-off value) G 0 is an "adjusted" parameter selected to fit the 

experimental data, and so as to prevent singularity in h when G-» 0 , and dev( ) is the deviatoric 
part of ( ). Note that eqs. (3.22) provide the simplest form of internal loading-unloading 
criterion; i.e., with the evolution for the internal strain variables, y tJ , during loading according 

to eq. (3.29b), but with a different dissipation function during unloading giving a “magnified” 
rate of evolution, y“ = to account for the “stiffen response upon 

unloading, with the moduli L u j]kl taken in the same form as for loading L tjkl in eq. (3.28a) with 
Go replacing G. However, altenative criterion are also possible to handle these observed 
readjustment of internal state [4], and may actually have direct impact on such cyclic-loading 
phenomena such as ratchetting. 


3.4 Generalized Viscoplasticity for Nonassociative Material 

As discussed in section 3.1, the nonassociative models are those which possess a 
partially (i.e., using only the fi-forms) or totally incomplete potential formulation. They 
either violate equations of state or evolution equations. Their popularity in the literature is 
mainly due to the added modeling flexibility to account for experimental phenomena by 
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using arbitrary rate forms in the flow and evolution equations. In this section, Freed’s 
viscoplastic model [24] is used as an example to discuss the theoretical basis of the 
nonassociative models. 

This viscoplastic model with creep and plasticity bounds was developed to 
represent the behavior of polycrystalline metals at high temperature, and is capable of 
accurately predicting short-term plastic strains, long-term creep strains and the 
interactions between them. It reduces analytically to a creep theory under steady-state 
conditions and becomes a plasticity theory at its rate-independent bound. This model is 
characterized by steady-state creep data, saturated hysteresis loops, and monatomic 
stress/strain curves. A general mathematical structure for viscoplasticity admits two 
tensorial internal variables, which are a short (a* ) and a long (a* ) back stress, and one 
scalar internal variable D, which is drag stress. In this model, the flow and evolution 
equations are no longer derived from Q-forms. Due to the nonassociative behavior of this 
model, it will be verified that the consistent viscoplastic-moduli matrix is no longer 
symmetric when the implementation of this model is discussed in Part II of this report. 

As mentioned, the discussion is limited to the case of an isotropic material under 
isothermal conditions. Freed’s original model is temperature dependent, and some material 
constants are temperature dependent. Under the present isothermal conditions, 
temperature effect is canceled from the original model equations, and those material 
constants which are the function of temperature are fixed during calculation at a given 
temperature. Tests may be given at different temperatures as long as the given temperature 
does not change during a test. In order to utilize the general framework of the algorithm 
represented in the previous section, the original tensor equations [78], which have the 
plasticity bound and influence of temperature variation are restricted and recast into the 
format used in the algorithm of GVIPS to facilitate the implementation of NAV. To be 
clear, the original equations of Freed’s model will be transformed into the new matrix 
format, (further details on these specific forms are given in Part II). 
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3.4.1. Elastic Behavior 

The stress, ay, is taken to be related to the infinitesimal total strain, ey, through the 
constitutive equations of an isotropic Hookean material, viz. 

Sy = 2/i(E ij ~£y ) where £& = 0 (3.33a) 

and 

0^=3^^ (3.33b) 

which are characterized by the shear, p, and bulk, k, elastic moduli, and where 

Sy = CJy - ^o-^j and Eij = e u -\ s u 5 ij (3.33c) 

denote the deviatoric stress and strain, respectively. Eq. (3.33a) characterizes the 
deviatoric stress response, while eq. (3.33b) characterizes the hydrostatic stress response. 

Based on the previous framework of section 3.3, the above equations of the elastic 
behavior are recast as, 


— Cijkl £ ld 

(3.34a) 

£ij = K+ e ij 

(3.34b) 


where are the same definition as before. The shear modules G and elastic modules 
E are defined as follows: 

G=|io + |iiT (3.35a) 

E = 2(l+v)G (3.35b) 

where T is temperature and v is Poisson’s ratio, po and pi are material constants. 
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3.4.2. Viscoplastic Behavior 


Plastic flow equation: 


’ n v = 


1 s —a 

— 1 v g 


2 s-a 


whose kinetics are governed by 


\\e p \=3Z if Its- all <K 


with the Von Mises norm ||s - a|| = J~( s o ~ ~ a u) 


(3.36) 


(3.37) 

(3.38) 


and where K > 0 denotes the plastic yield stress. 

The norms ( or magnitudes ) pertaining to the deviatoric tensors of this section are defined 
by 

n = ^ and (3 39) 

where Iy is any deviatoric “strain-like” tensor, and n ;j is any deviatoric “stress-like” 
tensor. 

Associated material functions quantifying plastic strain rate are as follows: 

when T t <T <T m 


exp 



when 0<T<T t 



(3.40) 
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Z = ^4sinh" 



(3-41) 


D 


Eq.3 .36 may be recast as follows: 


i’=f(J,D)T, ■ T,=M m ( a .-a H ) 


(3.42) 


where, f is the function of J and drag stress D, which is: 



(3.43) 


where J is quadratic function in terms of the invariants of the deviatoric components of the 


in which M.jti is defined in eq.(3.22). Due to the isotropic behavior assumed in the Freed’s 
original model, ^ = C, = 0, and Myki reduces to Pyki as defined in eq.(3.24a). All presently- 
used forms of this latter model are restricted to isotropic behavior; this is also the case 
here for all the tests performed later in section 4. 

The ( deviatoric, tensor-valued ) back stress, a, accounts for kinematic (flow- 
induced, anisotropic ) hardening effects. Here a is taken to be a finite sum of the 
individual back stresses, i.e., ay = a* + a' , where a* is called the short-range back 
stress, and a' ; is called the long-range back stress [24], The evolution of these back 
stresses is described through the deviatoric constitutive equations as follows: 




(3.44) 


Evolution equation: 


where H,, Hi are the short and long term hardening moduli, and L», Li are the thermal 
limiting states of the short and back stress, respectively, 


H, = |i and L, = f 


(C-D)(£>-A) 


(3.46a) 


SC 


H, = — and L, = (1-f) 


(C-D^D-D^) 


(3.46b) 


SC 


In eqs. (3.45), the first terms represent hardening while the second terms represent a 
deviation from strain hardening-a phenomenon called recovery. The evolution of the back 
stress accounts for the rapid change in stiffness that is observed during the transition from 
elastic to plastic behavior. 

Evolution equations (3.45a,b) for the back stress may be recast as follows: 



(3.47a) 



(3.47b) 


where, g s , g, are defined as: 



(3.48a) 
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g.= 


-^-«9ASinh n 

2L, 


rvn 


(3.48b) 


Another internal variable D, the drag strength, accounts for isotropic hardening 
effects. Because of its scalar nature, the strain hardening and dynamic recovery terms are 
combined into a single, nonlinear, hardening function, and the evolution of the drag 
strength is as follows: 


i>=h{p'\-Sr) ; ( 3 . 49 ) 

in which the first term of eq. (3.44) represents hardening, the second term recovery. The 
drag strength is bound by the interval D 0 <D< 

The hardening modules and other material constants for the drag strength are given as 
follows. 


h = h r 


(D-D 0 )/SC N 
sinh[(£>-Z> 0 )/<5C] 


y = Asinh’ 


' D-D 0 

SC 


Dnux = ^ ( C+D 0 ) 


h D = 


/», when T t <T<T m 

^ _VA r when 0<T<T t 

A 


(3.50) 


(3.51) 


(3.52) 


(3.53) 


This viscoplastic model is constrained such that creep is its lower ( steady-state ) 
bound, and plasticity is its upper ( rate-independent ) bound. When at steady state (i.e., 
a= 0, D= 0), this model reduces analytically to Odqvist’s creep theory, which greatly 
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simplifies the model and the process of its characterization. When below the plastic yield 
surface, i.e. F < 0, ( where, F = ||s-a|| - K or whenever ||s-a|| < K ), this model is a 
viscoplastic formulation. 

From the above equations it is noted that there are 13 material constants required 
to characterize this viscoplastic model: 

T m , the melting temperature; T t = TJ2 . 

Q , the activation energy for self-diffusion. 

& , the thermal function. 

H, , Hi are the hardening modules of the short, long-range back stress, respectively. 

f , the partitioning of back stress between short and long-range contributions. 

A, C, Do, 5, ho, hi,, m, n are other material constants used in the above eqs. 


3.5 The Implicit Integration Scheme 

Based on the theory discussed in the previous sections, the general implicit 
integration scheme is presented in this section. The details of the implementation and 
algorithm are discussed in Part II [38], As mentioned, the Euler scheme is presently used 
as the integrator, since it shows superior stability and convergence properties and has 
become the proven standard for the numerical integration of the viscoplastic rate 
equations. In this section, all equations and expressions are shown using a concise matrix 
notation for the integrated stress and internal stress fields. The contracted (Voigt) 
representations in vector forms for the components of the corresponding symmetric, 
second-order, tensor are utilized, with a slight abuse in notation to indicate the vector- 
matrix representations, that is, matrices are defined by under-curved symbols and vectors 
by underlined symbols. 
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3.5.1 General Form of Newton Iterative Scheme 


Consider a typical time step t„ -» tn+i, where the state variables at step n are 
known. Assume that At and Ae are given in the total strain field. The basic problem of 
evolution, then, is to update the state variables in a manner consistent with the governing 
equations; i.e., evolution equations and flow equations for the state variables. 


Using the implicit backward-Euler scheme, the update formulas are based on the 
following equation 


y. . — £ -f - d e * 

_n+l —n M “«+l 


(3.54) 


where is state variables, and dS n+I is the increment of the state variables, n is the step 
counter. For the GVIPS model, 


£.= 




(3.55a) 


whereas for the NAY model, 




yDj 


and the expression for the state variables increment is 


d = 


~ y 


— n+l 


(3.55b) 


(3.56) 


where AT Z is the iterative Jacobi matrix of state variables and R n+l is the residual function 
of state variables. The specific forms of and R n+} for the GVTPS and NAY models 


are defined in Part n. 
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To determine the updated values for the state variables, a local iterative solution is 
needed based on the eqs. (3.54, 3.56). For a typical iteration k -> k+1, 

?::: (3.57) 

note that d'L^ are evaluated based on the last updated states which means the iteration 
scheme is implicit. This equation represents the classical Newton-Raphson iteration 
scheme. It is well known that the basic Newton scheme is of “local” character; i.e., 
exhibiting excellent convergence properties (asymptotically quadratic rate) only when the 
starting trial solution is near the converged value, but for the trial solution which is far 
away from the converged value, it becomes very difficult to converge. Therefore, more 
robust numerical schemes are required for complex nonlinear problems. 


3.5.2 Line Search Scheme 

As mentioned, for highly-nonlinear problems, although the implicit scheme 
described above is unconditionally stable, its successful application still requires the use of 
an advanced numerical technique such as a line search method to produce an effective, 
robust solution algorithm. Its purpose is to guide the solution towards convergence, and 
is accomplied by searching for a scalar multiplier that adjusts the amount of the iterative 
increment vector to be updated within each iteration based on the optimization theory. 

The fine search method was utilized successfully in the global finite element 
solution of practical problems such as concrete cracking [20], and there is well-document 
evidence that the use of the line search is also essential for robust performance of global 
Newton’s method in inviscid elastoplasticity [68], At the global level, the concept of the 
line search algorithm pertains to minimizing the potential energy, that is, the work done by 
the residual force due to the iterative displacement. Here, and in this same spirit, the line 
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search method is firstly applied at the local (material point) level to drive the residual 
constitutive functions to be zero based on the optimization theory. 

Considering the line search method, eq. (3.57) need to be modified. Thus, for a 
typical iteration k -> k+1, the iterative procedure with line search may be described by the 
following expression. 


l£l=lL+nd£ 


7H-1 


(3.58) 


where 77 is a scalar ( 1 > 77 > 0 ) which adjusts the step size to optimize the iterative solution. 


Here, the objective of the line search method is to minimize 


1 


based on the 


optimization theory. Here, indicates the dot product of two vectors. In eq. (3.58), R 
is the “vector-valued” residual function of the state variables. For the GVIPS model, this 
is written symbolically as 


R k = 


U*J 


and for the NAY model, 




R k = 




(3.59a) 


(3.59b) 


For unconstrained problems, the residual comes from the difference between the current 
iterative value and the previous converged value. The calculations for dZ depend on the 
residual function and its derivatives at the previous step. The specific forms of R for 
GVIPS and NAY are defined in Part II [38], 


Now consider. 
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s(ri) = R(tj) • dl k (3 .60) 

in which E k and AS k are fixed, and R and s are functions of rj. When r\ = 0, 

So = s(r|=0) = R(q=0) • dZ k = Ro • dZ k (3.61) 


Ro is the residual function of the state variables at the end of the previous iteration. 
According to optimization theory, the best (optimum) solution is s(rj) = 0, but 
numerically, this is not realistically possible and it is also inefficient to try and achieve this 
objective. In practical problems, a ‘slack’ line search is used, with the aim of making the 
magnitude of S(rj) small in comparison with that of S(r|=0), i.e. 


s(ri = 0) 


<K 


(3.62) 


where 3u is the ‘line-search tolerance’. Based on past research experience, a suitable value 
for Pi. is of the order of 0.8 [20], After the successful search for -p is completed, eq. (3.58) 
keeps updating state variables until they satisfy the convergence condition. More details 
of the line search method are discussed in Part II [38], 


The convergence criterion here for Cauchy stress is in the following form. 


I -k+l „k 
CT n+ 1 °^n + 1 


i + l 
n+l I 


+ 0 . 01 £, 


< tol 


(3.63) 


where, I*! is the Euclidean norm of vectors, and tol is tolerance. Other internal variables also 
have similar convergence criterion. 
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4. Applications 

Considering the theoretical basis presented in the previous sections, a number of 
numerical tests are selected to demonstrate the capabilities, as well as various aspects, of 
the implicit iterative algorithm; e.g., validation, accuracy, robustness and computational 
efficiency. For any numerical method, accuracy, efficiency and a balance between the two 
are always major concerns. Validation of the algorithm may be proven by convergence and 
accuracy, while efficiency may be shown by CPU time comparisons with an explicit 
method or the same implicit method but with a different step size, but whose solution has 
the same accuracy. 

In this section, the performance of the implicit algorithm for both GVIPS and NAV is 
tested. For NAV, the new line search technique outlined earlier in section 3 is used to 
improve the efficiency and also guarantee the convergence. The details of the particular 
line search scheme utilized here in Part II of this report. The subincrementing technique is 
not used in the implementation of the NAV model. In the GVIPS model, there are regions 
of discontinuities in the state (a,j, ajj) space to account for stress-reversal effects, and for 
time steps during which these regions of discontinuity are traversed, the line search 
technique is no longer as powerful, as compared to the implementation of the NAV model 
(note that the line search is strictly valid only if the same “unique” forms of the residual 
functions apply during the considered step). Consequently, a subincrementing scheme is 
still needed to refine the global step size sufficiently to capture the point of discontinuity. 
For all other steps which are in the continuous region, the line search method will be 
shown to play an important role in terms of efficiency and convergence. 

For all test cases, the material parameters used in GVIPS model are associated 
with W/Kanthal (a unidirectional mmc model material), whereas the material parameters 
used in the NAV model are for copper. Material constants for both classes of viscoplastic 
models are given in Appendix 1. Consequently, no direct comparison between the GVIPS 
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and NAV models can be made as the different test cases are performed with these two 
materials. 

The solution schemes developed above were implemented in conjunction with a 4- 
node shell mixed finite element, as described in [66, 79], All numerical simulations 
reported herein were obtained using this element. 


4.1. Validation Tests 

Here, various tests are used to validate the implicit integration algorithm for both 
the NAV and GVIPS models. These tests include tension, cyclic, creep, and relaxation 
tests performed with different load steps. Additionally, a nonproportional multiaxial 
loading test was examined using the NAV and GVIPS model. Note that an explicit, one- 
step forward Euler integrator was used to obtain a converged reference solution for 
comparison with the implicit integrator. For all tests of GVIPS, the load is transverse to 
fibers. Among all tests, some tests are designed for assessing numerical performance of 
both models, and others are taken according to literature [36], 

Note that in the following tables, GIT is the average global iteration number, SUB 
is the average subincrement number, and LIT is the average local iteration number. Also, 
note that all of the CPU times are normalized with respect to the smallest CPU time. Here, 
CPU is the total time of a complete finite element analysis. Note that there are 8 
integration points for one homogeneous element. For the implementation method in the 
following tables, the implicit method denotes the backward Euler scheme, whereas the 
explicit method denotes the forward Euler. Finally, for convenience, all validation tests 
discussed subsequently are based on a finite element mesh consisting of a single element 
representing a material point. 
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4.1.1. Tension Tests 

This monotonic hardening test is used to demonstrate the accuracy and efficiency 
of the proposed implicit integration algorithm. For both NAV and GVIPS, the uniaxial, 
strain-controlled tensile test as shown in Fig. 1 was simulated. Where, the maximum 
applied strain was 0.20, and the total strain rate was 7.6 x 10^/s and 3.0 x 10' 5 /s for NAV 
and GVIPS respectively. 

4. 1 . 1 (a) NAV Isotropic Model 

Table 1 shows the results for a number of different load steps. First, note that the 
explicit integration method failed to reach the maximum strain level. After 10,000 steps, 
the explicit solution reached a maximum of 5% strain. To try and reach the maximum 
strain level, it was necessary to keep refining the step size, but this was stopped when the 
number of load steps and the CPU time became prohibitively large. From figure 2, the 
results of three different load step sizes are shown. The maximum stresses are 148, 148.57 
and 148.82 MPa for the 73, 20, and 10 load steps, respectively, and the differences are 
less than 1%, thus demonstrating the accuracy of the integration algorithm. It also may be 
concluded that the 73 step test is the converged solution by comparing its results with the 
explicit method in the region up to 5% strain. The agreement with the explicit reference 
solution is also validation of the implementation of the implicit algorithm. With regards to 
CPU time, it is found that the ratios of the explicit method, in the region up to 5% strain, 
and 73 step implicit method, in the region up to 5% strain, is 200/2.5=80, and the ratio of 
73 and 10 step cases is 3.2, which demonstrates the efficiency of the scheme. As an aside, 
the fact that the explicit solution required 10,000 steps is characteristic of explicit methods 
in general. For example, recent work by Arya [7] required from 40,000 to 100,000 
iterations. 
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a strain s 



Figure 1 Tension Schematic 


MPa 
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Table 1 Tension results for NAV 


Copper , Temperature=121 C° , g=7.6 x 1CT 6 /s , ^^= 0.20 , Dt=o=0.57 Mpa 


method 

number of 
load steps 

CPU time 

GIT 

LIT 

explicit 

10,000 

>200 

1 

0 

implicit 

73 

3.2 

2 

7 

implicit 

20 

1.1 

3 

8 

implicit 

10 

1 ( 7 ? *) 

3 

9 



Figure 2 Tension results for NAV (121 C° ) 
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4.1.1 (b) GVIPS Anisotropic Model 

Table 2 illustrates the results for three different load step increments, utilizing the 
GVIPS model. Once again, 10,000 steps were applied for the explicit solution, here 
however, the explicit solution reached the maximum strain. The results of three different 
cases are the same as shown in fig. 3, thus demonstrating the accuracy of the integration 
algorithm and also the validation of the implementation of the implicit algorithm. The 
conclusion may be drawn that the 10 step case is the converged solution by comparing its 
results with the explicit method. With regards to CPU time, it is found that the ratios of 
the explicit method and the 10 step implicit method is 312, and the ratio of the 100 and 10 
step case is 2.0, thus demonstrating the efficiency of the scheme. 


4.1.2. Cyclic Loading Tests 

This example serves to assess the accuracy of the stress predictions and to 
demonstrate the robustness of the implicit algorithm. For both NAV and GVIPS, the 
uniaxial, cyclic (non-monotonic), strain-controlled test, shown in Fig. 4, was simulated. 
The maximum strain amplitude was ±0. 144 percent applied at a rate of 0.002/s. 

4. 1 .2(a) NAV Isotropic Model 

As shown in figure 5, the 100, 10, and 5 step cases agree reasonably with the 
converged explicit solution. For the 5 step case, the stress is updated from approximately - 
220 MPa to 200 MPa within one step, which demonstrates the robustness of the 
algorithm. Comparing the 100 step case with the explicit method which needs 10,000 
global steps for one complete cycle, the CPU time ratio is approximately 340/6^60 (table 
3). Considering the comparable accuracy of these two cases, a CPU ratio of 60 
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Table 2 Tension results for GVIPS 


W/Kanthal , g=3.0x 10' 5 /s , £^=0.20 


method 

number of 
load steps 

CPU time 

GIT 

LIT 

explicit 

10,000 

312 

1 

0 

implicit 

100 

2.0 

2 

3 

implicit 

10 

1.0 (73 s) 

3 

4 



Strain 


Figure 3 Tension results for GVIPS 
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Figure 4 Cyclic Schematic 
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Table 3 Cyclic results for NAY 


Copper , Temperature=21 C° , g =2x1 O' 3 /s , £•^=0.0144 , Dt=o=5.9 MPa 


method 

number of 
load steps 

CPU time 

GIT 

LIT 

explicit 

10,000 

340 

3 

0 

implicit 

100 

6.0 

2 

4 

implicit 

10 

1.05 

3 

10 

implicit 

5 

1 ( 38 s) 

3 

22 



Strain 


Figure 5 Cyclic results for NAV ( 21 C° ) 
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demonstrates a significant improvement in efficiency and may be attributed to the use of 
the line search method. 

4. 1 .2(b) GVIPS Anisotropic Model 

Figure 6 shows the explicit converged solution as compared to the implicit 100, 10 
and 5 step cases. The fact that all of the implicit cases match the converged solution 
demonstrates the accuracy, efficiency and robustness of the algorithm. In addition note 
that even the 5 step case matches the converged solution at all points of the hysteresis 
loop. The previously discussed discontinuities within the GVIPS formulation occurs along 
the path from A to B and from C to D as shown in figure 4. Those steps which contain the 
discontinuity point require subincrementing in order to cross the discontinuity. Taking the 
10 step case as an example, the third and seventh global steps contain the discontinuity. 
Thus, only these two steps required subincrements which averaged 22 in this case. As 
such, comparing the 10 step case with the explicit solution (table 4), the CPU ratio is 
approximately 1 80 which demonstrates a marked saving in CPU time yet with comparable 
accuracy. 


4.1.3. Creep Tests 

The creep test is used to test both the accuracy and robustness (stability) of the 
implicit integration algorithm. In a creep test, the primary creep region is controlled by 
accuracy due to the high inelastic strain rates involved. In the secondary (steady state) 
creep region, the stability or robustness of the integration algorithm is tested. Both the 
NAV and GVIPS models will be used in the creep tests. For the NAV model one creep 
test of 1000 seconds will be performed, at a constant creep stress of 24 Mpa, whereas for 
the GVIPS model, one creep test of 1800 seconds at a creep stress of 70 MPa will be 
performed. Note in Fig. 7, e 0 is the total strain at the beginning of the hold time and s f is 
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Table 4 Cyclic results for GVIPS 


W/Kanthal , g=2xl0' 3 /s , £^=0.0144 


method 

number of 
load steps 

CPU time 

GET 

LIT 

explicit 

10,000 

180.0 

3 

0 

implicit 

100 

5.0 

2 

4 

implicit 

10 

1.05 

4 

10 

implicit 

5 

1 ( 54 s) 

10 

20 



Figure 6 Cyclic results for GVIPS 
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Figure 7 Creep Schematic 
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the amount of total strain at the end of the hold time. A large initial loading rate was 
chosen for both creep tests such that the initial load-up region may be considered as a 
purely elastic loading without any inelastic strain. 

4. 1.3(a) NAV Isotropic Model 

The explicit method with a total of 20,400 global steps was determine to have 
converged and is used as the reference solution (Fig. 8). The implicit method, was run 
using 480 and 880 steps (table 5). With the 880 step case producing a solution that differs 
from the converged solution by approximately 1%. Note that it is in the primary creep 
region where accuracy is most important. It is well established that a larger number of 
steps are required to satisfy the requirement of high accuracy for the rate-dependent 
problem. With regards to CPU time, upon comparing the explicit method with the 880 
step implicit case, there is a ratio of 32/1.6=20 which is a significant reduction in 
computation time. 

4.1.3(b) GVIPS Anisotropic Model 

The converged solution was obtained using the explicit method with a total of 
5840 global steps (Fig. 9). Note that the difference between the explicit method and 805 
step implicit method is less than 2% (table 6). Comparing CPU times, the ratio between 
the explicit method and the 805 step case is approximately 6.0, which demonstrates the 
efficiency of the implicit method. The 405 step case has a comparatively large error in 
comparison with the explicit solution. But, note that the secondary creep rate of the 405 
step case does appear to match that of the converged solution. As expected, since the 
error may be attributed to the step size (accuracy and not stability), an improved primary 
creep response was predicted (805 step case) thus basically shifting the curve upward and 
closer to the converged solution. 
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Table 5 Creep results for NAY 


method 

7 w — r — 

number of 
load steps 

CPU 

«. (10°) 

- — 
«/ (10- 5 ) 

GIT 

, — L-U ~ 

LIT 

creep time 

(s.) 

explicit 

20,400 

32 

0.6273 

0.8114 

1 

0 

1000 

implicit 

880 

1.6 

0.6313 

0.8022 

2 

1 

1000 

implicit 

480 

1 (954s) 

0.6313 

0.7853 

3 

2 

1000 



Creep time (hours) 


Figure 8 Creep results for NAV (121 C° ) 
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Table 6 Results of creep tests for GVIPS 


W/Kanthal , cq = 10 


csi, cr = 1 ksi/s 


method 

number of 
load steps 

CPU 

(s.) 

(10* 3 ) 

«, (10- 3 ) 

GIT 

LIT 

creep time 
(sec.) 

implicit 

405 

1 (700) 

0.3850 

1.1016 

2 

2 

1800 

implicit 

805 

1.7 

0.3850 

1.2630 

2 

1 

1800 

explicit 

5840 

11.0 

0.3759 

1.2884 

2 

0 

1800 



Figure 9 Creep results for GVIPS 
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4.1.4. Relaxation Tests 

The stress relaxation test may also be used to assess the accuracy of the integration 
method. Subsequent to load up, typically, the relaxation test. Fig. 10 is characterized by a 
marked decrease in the initial stress at a high rate. The amount of stress relaxation and 
corresponding rate is dependent on the specific material model. Thus, an accurate 
integration method is required so as to properly predict the highly non-linear stress 
response. Relaxation tests were performed for both the NAV and GVIPS models. For the 
NAV model, the constant strain is 0.01 and the relaxation time for the test is 1000 
seconds. Whereas for GVIPS, the constant strain is also 0.01, but the relaxation time is 
1800 seconds. 

4. 1 .4(a) NAV Isotropic Model 

The converged, reference solution for the 1,000 second relaxation test using the 
explicit method required a total of 80,000 steps (table 7). This large number of steps was 
required to accurately predict the highly non-linear stress response. The implicit method 
was run using 110 and 210 steps, with the 210 step case essentially matching the 
converged solution (Fig. 1 1). Again, note that it is in the region where the stress decreases 
at a high rate that the solution is most sensitive to accuracy. Comparing the CPU time 
between the explicit and the 210 step implicit, there is a ratio of 900 in CPU time. Thus 
for comparable accuracy, the implicit method is far more efficient. As an aside, note that 
there is no significant predicted stress relaxation (i.e. < 4 MPa) for Freed’s model. This 
may be due to the simplifying assumption that there is no explicit thermal recovery of the 
short- and long-range back stresses in his model. 

4.1.4(b) GVIPS Anisotropic Model 

The converged solution using the explicit method required a total of 3894 global 
steps (table 8). Fig. 12 shows the results for 22, 52 and 102 steps. Note that it is in the 
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Figure 10 Relaxation Schematic 


Table 7 Relaxation results for NAY 
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Copper , Temperature=3 16 C° , e f = 0.01 , £=lxl0' 3 s' 1 , Dt-o=0.50MPa 


method 

number of 
load steps 

CPU 

O’o 

MPa 

°7 

MPa 

GIT 

LIT 

relaxation 
time (s) 

explicit 

80,000 

1350 

36.331 

32.848 

1 

0 

1000 

implicit 

210 

1.5 

36.248 

32.879 

2 

1 

1000 

implicit 

110 

1 (160s) 

36.248 

33.093 

2 

2 

1000 



Relaxation time (hours) 


Figure 1 1 Relaxation results for NAV (3 16 C° ) 
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Table 8 Results of relaxation tests for GVIPS 


W/Kanthal ,e f = 0.01 , e =lxl0' 3 s' 1 


method 

number of 
load steps 

CPU 

(s) 

1 

o 

b 

<j f ksi 

GIT 

LIT 

relaxation 
time (s) 

implicit 

52 

1(71) 

49.887 

24.316 

2 

2 

1800 

implicit 

102 

1.8 

49.887 

23.625 

2 

1 

1800 

explicit 

3894 

53.9 

50.701 

23.430 

2 

0 

1800 



Figure 12 Relaxation results for GVIPS 
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initial relaxation region, (time < 0.2 hrs.) in which the solution is most sensitive to time 
step size and thus accuracy. The difference between the explicit method and the 102 step 
implicit case is less than 1% at the end of the relaxation period and approximately 8% 
difference at 0.05 hrs, where the solution is highly non-linear. Again, comparing the CPU 
time between the explicit and 102 step implicit methods, it is shown, see table 8, that a 
speed up ratio of approximate 30 times is indicated. Once again illustrating the efficiency 
of the implicit method. Finally, note that GVIPS predicts a significant amount of stress 
relaxation as suggested by experimental evidence, that is, approximately a 50% reduction 
in the initial stress after load up. 

4.1.5 Nonproportional Multiaxial Tests 

This final test examines the integration algorithm under nonproportional loading 
conditions. The element is subjected to the nonproportional strain path shown on Fig. 13. 
This path is taken from the experimental work of Lamba and Sidebottom [36], This 
applied strain path was chosen because of its wide range of angles between segments with 
equal positive and negative peaks. The 8 paths combine shear, axial loading and 
unloading. The end points of the path segments were numbered from 0 to 8 and the path 
was traversed in a numerically increasing sequence. Note that the experiment given in Ref. 
[36] was conducted on a saturated specimen, that is, complete stabilization with regards to 
hardening has occurred. 

4. 1 .5(a) NAV Isotropic Model 

The stress response of the model to the applied nonproportional strain path 
appears in Fig. 14. This test demonstrates the robustness of the scheme discussed. The 
predicted response, shown in Figure 14, matches the experimental results (not shown in 
fig. 14) given in Ref. [36] very well. In addition, note that even one step per path ( figure 
14 ) can successfully complete this complicated loading history. Because the 50 step and 
10 step cases give identical results, it is concluded that 50 steps has converged to the 
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Figure 13 Nonproportional loading path 
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Figure 14 Nonproportional loading results for NAV (121 C° ) 
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correct solution. Finally the CPU time for 1, 10, 50 steps per path are 1 ( 35s ), 3.8 and 
21.4, respectively. 

4.1.5(b) GVIPS Isotropic Model 

This GVIPS model was originally used to characterize W/Kanthal whose material 
constants are given in Appendix 1. In order to compare the stress response of the GVIPS 
model under nonproportional strain path with the experiment given in [36], the material 
constants needed to be changed to fit the copper material used in the experiment. The 
new material constants are given in Appendix 2, wherein the material is taken to be 
isotropic, i.e., co=r|=1.0. Based on the new material constants, the predicted response of 
the GVIPS model is shown in Fig. 15, with the 10 step and 100 step cases giving almost 
identical results. Given the very crude characterization of the GVIPS model (essentially 
fitting the limit state in a single pure shear test), it is remarkable how accurate it represents 
the experimental results. 


4.2. Efficiency : Comparative Studies 

In addition to validating the proposed integration method, efficiency studies are 
necessary since efficiency is a very important feature that is required in order to solve 
large, realistic finite element problems. For such practical problems, it is desirable and 
necessary to obtain an accurate solution in the most economical manner, that is, using the 
minimum CPU time. 

Many mathematical techniques have been applied in order to develop computationally 
efficient methodologies for nonlinear analyses in which a large global time step size may 
be chosen in order to reach convergence with sufficient accuracy. Subincrementing has 
been shown to be successful in guaranteeing convergence when using a large global time 
step. The basic idea in subincrementing is to refine the global step size to achieve local 
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Figure 15 Nonproportional loading results for GVIPS 
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iteration convergence, as well as maintaining sufficient accuracy with regards to the 
updated fields, i.e. internal, state variables, stress, strain. For difficult load histories such as 
the cyclic test or for a complex, numerically stiff constitutive model, such as Freed’s NAV, 
too many subincrements may be required, which directly results in large CPU times. 
Basically, subincrementing is useful for improved convergence, but it is not efficient since 
it does not improve the rate of convergence. Also, to date, there seems to be no 
agreement with regard to a single criterion of general applicability for determing the 
“optimal” number of subincrements. With these points in mind, a line search method was 
introduced here as a viable, and far more effective, alternative in providing improved 
convergence and solution economy. A detailed discussion of the line search method will be 
presented in Part II [38], 

In order to compare the efficiency between the line search method and the 
subincrementing method, all the tests described below were performed using both of the 
two methods. In the following. Method 1 utilizes subincrementing only, while Method 2 
u tiliz es only the line search method. Recall that for GVIPS, there are regions of 
discontinuities in the state space to account for stress-reversal effects which are loading 
and unloading conditions. In these steps, subincrementing is combined with line search in 
order to refine the step size to locate the point of discontinuity. This is required since the 
premise of line search is that equations optimized are continuous. Thus, Method 2 consists 
of a combination of subincrementing and line search methods, when the GVIPS model is 
utilized under cyclic loading conditions. For either the NAV or GVEPS model under 
monotonic loading condition. Method 2 is a line search only method. The results 
presented in tables 9 and 10 demonstrate the advantage of using a line search in the 
integration algorithm. 

Comparing method 1 and 2 in the above tables, note that the CPU time of method 
1 is always larger than the CPU time of method 2. In particular for GVIPS, the CPU ratios 
are 6.0, 4.5, 1.5, and 1.1 for the tension, cyclic, creep, and relaxation tests, respectively. 
For NAV, the CPU ratios are 13.65, 14.5, 1.55, and 22.2 for the tension, cyclic, creep. 
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Table 9 Results of comparision for GVIPS 
Table 9a tension test, global steps=10 £,^=0.20 


Method 

CPU ( seconds) 

GIT 

SUB 

LIT 

subincrementing 

6.0 

11 

11 

3 

line search * 

1 (73) 

3 

0 

4 


Table 9b cyclic test, global steps=25 €mix=0.0144 


Method 

CPU ( seconds) 

GIT 

SUB 

LIT 

subincrementing 

4.5 

3 

8 

2 

line search * 

1 (33.91) 

2 

0 

4 


Table 9c creep test, global steps= 105 o=70MPa 


Method 

CPU ( seconds) 

GIT 

SUB 

LIT 

subincrementing 

1.5 

4 

4 

2 

line search * 

1 ( 178 ) 

4 

0 

2 


Table 9d relaxation test, global steps=105 €nux = 0.01 


Method 

CPU ( seconds) 

GIT 

SUB 

LIT 

subincrementing 

1.1 

3 

2 

2 

line search * 

1 ( 122 ) 

3 

0 

2 


* Combination of the line search method and subincrementing (only when crossing 
regions of discontinuity). 
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Table 10 Results of comparision for NAV 

Table 10a tension test, global steps= 10 £^=0-20 


Method 

CPU ( seconds) 

GIT 

SUB 

LIT 

subincrementing 

13.65 

16 

15 

3 

line search 

1 (77) 

3 

0 

9 


Table 10b cyclic test, global steps=5 £max = 0.0144 


Method 

CPU ( seconds) 

GIT 

SUB 

LIT 

subincrementing 

14.5 

8 

15 

5 

line search 

1 (38) 

3 

0 

22 


Table 10c creep test global steps=23 o= 70 Mpa 


Method 

CPU ( seconds) 

GIT 

SUB 

LIT 

subincrementing 

1.55 

25 

15 

5 

line search 

1 (768 ) 

20 

0 

15 


Table lOd relaxation test, global steps=22 £max = 0.01 


Method 

CPU ( seconds) 

GIT 

SUB 

LIT 

subincrementing 

22.2 

6 

24 

2 

line search 

1 (50) 

2 

0 

3 
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and relaxation tests, respectively. For NAV the largest CPU ratio is 22.2, i.e., method 2 
with line search is approximately 95% faster than method 1 with subincrementing. This 
shows a significant improvement in efficiency particularly for a problem which requires 
accuracy, as in the case of the relaxation test. With regards to the number of iterations 
both globally and locally, line search usually required more local iterations, but this is due 
to the fact that it does not use any subincrements. From the results, it is apparent that the 
small additional local calculations required for the line search method leads to a saving in 
the total calculations. Specifically, some tests may need a large number of subincrements. 
For example, in the relaxation test of NAV, the number of subincrements for method 1 is 
24 with 2 local iterations, table 10b. On the other hand, with line search, method 2, only 3 
local iterations are required, thus resulting in a significant reduction in CPU time. 

Finally, comparing tables 9, and 10, notice that for GVIPS, the CPU ratios of 
method 1 to method 2 are always less than NAV’s. This trend may indicate that it is more 
difficult to integrate NAV without the help of sophisticated numerical techniques. On the 
other hand, since line search does not appear to significantly affect local convergence, 
GVIPS may be comparably easier to integrate and may only require a simple 
subincrementing scheme. 

4.3. A General Structural Problem Using a NAV Isotropic Model 

As has been stated throughout this report, the ultimate test of the integration 
method is to demonstrate its accuracy, robustness and efficiency when applied to a 
realistic structural analysis problem. In this context, a perforated plate as presented in 
reference [29] and [76], was analyzed using Freed’s NAV model. A quarter of the plate is 
discretized for FE analysis (see figure 16). A linearly increasing specified displacement, 
with a maximum value of 0.5 percent ( Al/w ) at a rate of 6.67xl0' 5 /s, was applied along 
the upper boundary of the plate. Material constants are the same as those given in 
Appendix 1. The results shown in Figure 17-19 are the effective stress distribution at the 


Figure 16 Structure mesh 
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final load level when 1, 5 and 20 global steps are utilized respectively. Also in order to 
check the convergence and efficiency of the implicit integration scheme, an explicit test is 
run by using 5000 global steps as the reference solution. The result shown in Figure 20 
and Figure 21 are the effective inelastic strain distribution at the final load level by using 
the implicit scheme ( 20 global steps ) and explicit scheme ( 5000 global steps ), 
respectively. Note that the difference between these two figures is less than 1% which 
means the 20 step case converged and may be chosen as the reference solution. 
Comparing the 5 and 20 step cases (Figure 18 and 19), note that the effective stress 
distribution of the 5 step case is qualitatively close to the 20 step case which shows 
reasonable accuracy for even the 5 step case. Also, convergence was achieved using only 
one global step, this demonstrates the robustness of the implicit integration algorithm. 
Although the accuracy has decreased (12% difference in maximum J 2 ) which is 
understandable since this is a first order integration method and just one step is used. The 
CPU time for the explicit scheme, and the 20, 5, and 1 step implicit cases are 12570, 251, 
91 and 25 minutes, respectively. The CPU ratio of the explicit versus the implicit scheme ( 
20 steps ) is more than 50, thus showing the significant efficiency of the implicit 
integration scheme. Comparing the implicit scheme itself, the CPU ratio of the 20 step 
case to that of the 5 step case is approximately 2.7. This a significant saving in 
computation time yet the results are almost identical, i.e. the difference is less than 2%. 
Thus, for the complex structural problem, convergence may be reached in an efficient 
manner without a significant loss in accuracy. Again, there will always be a balance 
between solution accuracy and computational efficiency. 


5. Summary & Conclusions 


The general computational framework described in this report using the backward 
Euler fully-implicit integration method has been determined to be successful for both the 


Figure 17 Effective stress distribution (1 step) 
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Generalized Viscoplasticity with Potential Structure (GVIPS) and the Non-Associative 
Viscoplastic models (NAV). Only these Newton iterative schemes will provide, uniformly- 
valid, convergent, implicit integrations for both GVIPS and NAV classes of viscoplastic 
models. The algorithm was written in a concise matrix format so as to provide sufficient 
generality and is automatically valid for both isotropic and anisotropic cases in full space 
or subspaces. 

All numerical tests show that the use of the implicit integration method is robust, 
unconditionally-stable and provides first-order accuracy for both “small” and “large” At 
that is, the integration method provides a balanced combination of stability and accuracy. 
Even with nonproportional, multiaxial stress-strain states, sufficient accuracy and good 
convergence characteristics were demonstrated for the implemented scheme. 

In addition, it has been generally observed that explicit schemes either fail 
completely or become prohibitively expensive, especially for the NAV class of models. 
Specifically the time step size must be inherently very small in order to guarantee 
convergence, which is not efficient. It is believed that this point has been clearly 
demonstrated in the results presented here for the explicit integration cases. That is, it was 
required that the number of time steps be on the order of 10,000 to 80,000 in order to 
obtain a solution. These results appear to be consistent with other recent work, for 
example that of Arya [7], that uses various explicit schemes. 

The implicit integration scheme with simple subincrementing, which improves 
convergence, has been shown to be efficient when used for GVIPS and NAV. But on the 
other hand, subincrementing alone is not efficient. That is, many subincrements and local 
iterations are required. As a result, the ’’slack” method, which is the form of line search in 
the present algorithm, enables the NAV class of models to converge stably with sufficient 
accuracy and significantly improved efficiency. 
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In comparing the recently-developed GVIPS class model versus some forms of the 
presently-used NAV models found in the literature, the GVIPS formulation posses a more 
theoretically sound basis from the thermodynamics standpoint, and exhibit symmetries 
which are desirable with respect to computation in the corresponding integrated forms of 
their flow and evolution equations. For example, for global iterations, GVIPS provide a 
fully-consistent tangent stiffness, but for NAV, only an approximate tangent stiffness is 
obtained. From the numerical tests in this report, it was found that NAV exhibited greater 
difficulty in integration when compared to the GVIPS model. That is, NAV required more 
iterations in order to reach a solution. However, by means of a line search scheme, NAV 
(because of its continuous formulation), can now be integrated more easily even without 
subincrementing. The overall efficiency of integrating the GVTPS class of viscoplastic 
models is even more impressive when one recalls that a combination of subincrementing 
and line search methods must be applied; because of the discontinuity in GVIPS evolution 
equations, due to the need to correctly account for unloading. In particular, in the 
discontinuous regions, subincrementing is still required to refine the time steps sufficiently 
in order to capture the discontinuity, but is only needed in those areas (few steps) near the 
discontinuity. Alternatively, in the continuous regions, of GVIPS the use of a line search 
technique applies and does contribute significantly to improve convergence. Thus, 
intelligent combination of line search and subincrementing guidance convergence as well 
as overall efficiency for the GVIPS class of models. 
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Appendix 1 : Material constants 


NAY material constants ( Copper ): 


Constant 

Value 

Units 

T t 

678 

K 

Mo 

43000 

MPa 

Mi 

-17 

MPa/K 

ho 

50 

MPa 

hi 

15 

MPa 

Q 

200,000 

J/mol. 

K 

8.314 

J/mol. K 

A 

20,000,000 

1/s 

n 

4.5 


c 

13 

MPa 

Do 

0.13 

MPa 

f 

0.75 


6 

0.035 


m 

0.5 


Ho 

20 


V 

0.36 



GVIPS material constants (W/Kanthal): 


Constant 

Value 

Units 

n 

1.5 


M 

2.5E4 

hr 

m 

1.5 



0.5 


R 

7.0E-5 

MPa/hr 

H 

12.6E4 

MPa 

Go 

0.05 


CD 

2.7 


_Q_ 

1.0 


kt 

5.6 

MPa 
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Appendix 2: Modified material constants 


GVIPS material constants for nonproportional test: 


Constant 

Value 

Units 

n 

2.5 


_LL_ 

8.1238E12 

hr 

m 

3.0 


j_ 

2.5 


R 

7.0E-4 

MPa/hr 

H 

2.0E4 

MPa 

Go 

2.5 


CO 

1.0 


n 

1.0 


kt 

0.6823 

MPa 
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